home *** CD-ROM | disk | FTP | other *** search
/ United Public Domain Gold 2 / United Public Domain Gold 2.iso / utilities / pu309.dms / pu309.adf / Sat_Tracking / sat-partb.c < prev    next >
C/C++ Source or Header  |  1990-09-29  |  7KB  |  226 lines
  1. /* N3EMO Orbit Simulator routines
  2.    Copyright (C) 1986 Robert W. Berger  */
  3.  
  4. #include "df0:include/stdio.h"
  5. #include "df0:include/math.h"
  6. #define PI 3.14159265
  7. #define PI2 (PI*2)
  8. #define MinutesPerDay (24*60.0)
  9. #define SecondsPerDay (60*MinutesPerDay)
  10. #define HalfSecond (0.5/SecondsPerDay)
  11. #define EarthRadius 6378.16      /* Kilometers */
  12. #define EarthFlat (1/298.25)      /* Earth Flattening Coeff. */
  13. #define SiderealSolar 1.0027379093
  14. #define DegreesPerRadian (180/PI)
  15. #define RadiansPerDegree (PI/180)
  16. #define ABS(x) ((x) < 0 ? (-(x)) : (x))
  17. #define SQR(x) ((x)*(x))
  18. #define RefYear 79
  19. #define RefSidereal 0.27610467
  20. #define RefDay 0      /* Day of week */
  21. char *MonthNames[12] = { "Jan","Feb","Mar","Apr","May","Jun","Jul",
  22.          "Aug","Sep","Oct","Nov","Dec" };
  23. int MonthDays[12] = {31,28,31,30,31,30,31,31,30,31,30,31};
  24. char *DayNames[7] = { "Sunday","Monday","Tuesday","Wednesday","Thursday",
  25.          "Friday","Saturday"};
  26. /* Solve Kepler's equation               */
  27. /* Inputs:                     */
  28. /*   MeanAnomaly   Time since last perigee, in radians.   */
  29. /*         PI2 = one complete orbit.      */
  30. /*   Eccentricity   Eccentricity of orbit's ellipse.   */
  31. /* Output:                     */
  32. /*   TrueAnomaly   Angle between perigee, geocenter, and   */
  33. /*         current position.         */
  34. double Kepler(MeanAnomaly,Eccentricity)
  35. register double MeanAnomaly,Eccentricity;
  36. {
  37. register double E;      /* Eccentric Anomaly         */
  38. register double Error;
  39. register double TrueAnomaly;
  40.     
  41.     E = MeanAnomaly;   /* Initial guess */
  42.     do
  43.         {
  44.    Error = (E - Eccentricity*sin(E) - MeanAnomaly) 
  45.       / (1 - Eccentricity*cos(E));
  46.    E -= Error;
  47.    }
  48.    while (ABS(Error) >= 0.000001);
  49.     
  50.     if (ABS(E-PI) < 0.000001)
  51.         TrueAnomaly = PI;
  52.       else
  53.         TrueAnomaly = 2*atan(sqrt((1+Eccentricity)/(1-Eccentricity))
  54.             *tan(E/2));
  55.     if (TrueAnomaly < 0)
  56.    TrueAnomaly += PI2;
  57.     return TrueAnomaly;
  58. }
  59. GetSubSatPoint(EpochTime,EpochRAAN,EpochArgPerigee,Inclination,
  60.    Eccentricity,RAANPrecession,PerigeePrecession,Time,TrueAnomaly,
  61.       Latitude,Longitude)
  62. double EpochTime,EpochRAAN,EpochArgPerigee,Inclination,Eccentricity;
  63. double RAANPrecession,PerigeePrecession,Time;
  64. double TrueAnomaly,*Longitude,*Latitude;
  65. {
  66.     double RAAN,ArgPerigee;
  67.     int i;
  68.     ArgPerigee = EpochArgPerigee + (Time-EpochTime)*PerigeePrecession;
  69.     *Latitude = asin(sin(Inclination)*sin(ArgPerigee+TrueAnomaly));
  70.     RAAN = EpochRAAN - (Time-EpochTime)*RAANPrecession;
  71.     *Longitude = -acos(cos(TrueAnomaly+ArgPerigee)/cos(*Latitude));
  72.     if ((Inclination > PI/2 && *Latitude < 0) 
  73.           || (Inclination < PI/2 && *Latitude > 0))
  74.    *Longitude = -*Longitude;
  75.     *Longitude += RAAN;
  76.     /* Convert from celestial to terrestrial coordinates */
  77.     *Longitude -= PI2*(Time*SiderealSolar + RefSidereal);
  78.     /* Make West be positive */
  79.     *Longitude = -*Longitude;
  80.     /* i = floor(Longitude/2*pi)   */
  81.     i = *Longitude/PI2;
  82.     if(i < 0)
  83.    i--;
  84.     *Longitude -= i*PI2;
  85. }
  86.  
  87. GetPrecession(SemiMajorAxis,Eccentricity,Inclination,
  88.    RAANPrecession,PerigeePrecession)
  89. double SemiMajorAxis,Eccentricity,Inclination;
  90. double *RAANPrecession,*PerigeePrecession;
  91. {
  92.   *RAANPrecession = 9.95*pow(EarthRadius/SemiMajorAxis,3.5) * cos(Inclination)
  93.        / SQR(1-SQR(Eccentricity)) * RadiansPerDegree;
  94.   *PerigeePrecession = 4.97*pow(EarthRadius/SemiMajorAxis,3.5)
  95.       * (5*SQR(cos(Inclination))-1)
  96.        / SQR(1-SQR(Eccentricity)) * RadiansPerDegree;
  97. }
  98. GetBearings(SiteLat,SiteLong,SiteAltitude,
  99.    SatLat,SatLong,Radius,Azimuth,Elevation,Range)
  100. double SiteLat,SiteLong,SiteAltitude,SatLat,SatLong,Radius;
  101. double *Azimuth,*Elevation,*Range;
  102. {
  103.     double SinSiteLat,sinSiteLong,SinSatLat,SinSatLong;
  104.     double CosSiteLat,cosSiteLong,CosSatLat,CosSatLong;
  105.     double CosBeta,SinBeta;
  106.     double LongDiff;
  107.     SinSiteLat = sin(SiteLat); sinSiteLong = sin(SiteLong);
  108.     CosSiteLat = cos(SiteLat); cosSiteLong = cos(SiteLong);    
  109.     SinSatLat = sin(SatLat); SinSatLong = sin(SatLong);
  110.     CosSatLat = cos(SatLat); CosSatLong = cos(SatLong);    
  111.     LongDiff = SiteLong - SatLong;
  112.     CosBeta = SinSiteLat*SinSatLat+CosSiteLat*CosSatLat*cos(LongDiff);
  113.     SinBeta = sqrt(1-SQR(CosBeta));
  114.     *Azimuth = acos((SinSatLat- SinSiteLat*CosBeta)/CosSiteLat/SinBeta);
  115.     if (LongDiff < -PI)
  116.         LongDiff += PI2;
  117.     if (LongDiff > PI)
  118.         LongDiff -= PI2;
  119.     if (LongDiff < 0)
  120.    *Azimuth = PI2 - *Azimuth;
  121.     /* Convert to geocentric radius */
  122.     SiteAltitude += EarthRadius*(1-EarthFlat/2+EarthFlat/2*cos(2*SiteLat));
  123.     *Elevation = atan((Radius*CosBeta-(SiteAltitude))
  124.              /(Radius*SinBeta));
  125.     *Range = sqrt(SQR(Radius) + SQR(SiteAltitude)
  126.               -2*Radius*SiteAltitude*CosBeta);
  127. }
  128.  
  129. PrintDate(OutFile,Day)
  130. FILE *OutFile;
  131. {
  132.     int Month,Year,DayOfWeek;
  133.     DayOfWeek = (Day-RefDay) % 7;
  134.     Year = RefYear;
  135.     while (Day <= 0)
  136.    {
  137.    Year -= 1;
  138.    if (Year%4 == 0)
  139.        Day += 366;
  140.     else
  141.        Day += 365;
  142.    }
  143.     while ((Year%4 == 0 && Day > 366) || (Year%4 != 0 && Day > 365))
  144.    {
  145.    if (Year%4 == 0)
  146.        Day -= 366;
  147.     else
  148.        Day -= 365;
  149.    Year += 1;
  150.    }
  151.     
  152.     if (Year%4 == 0)
  153.         MonthDays[1] += 1;   /* Leap year */
  154.     
  155.     Month = 0; 
  156.     while (Day > MonthDays[Month])
  157.         {
  158.    Day -= MonthDays[Month];
  159.    Month += 1;
  160.    }
  161.     fprintf(OutFile,"%s  %d %s %d",DayNames[DayOfWeek],Day,
  162.           MonthNames[Month],Year);
  163.     if (Year%4 == 0)
  164.         MonthDays[1] -= 1;   /* Leap year */
  165. }
  166.  
  167. PrintTime(OutFile,Time)
  168. FILE *OutFile;
  169. double Time;
  170. {
  171.     int day,hours,minutes,seconds;
  172.     day = Time;
  173.     Time -= day;
  174.     if (Time < 0)
  175.         Time += 1.0;   /* Correct for truncation problems with negatives */
  176.     hours = Time*24;
  177.     Time -=  hours/24.0;
  178.     
  179.     minutes = Time*MinutesPerDay;
  180.     Time -= minutes/MinutesPerDay;
  181.     
  182.     seconds = Time*SecondsPerDay;
  183.     seconds -= seconds/SecondsPerDay;
  184.           
  185.     fprintf(OutFile,"%02d%02d:%02d",hours,minutes,seconds);
  186. }    
  187. /* Get the Dayn(SatLong);
  188.     CosSatLat = cos(SatLat); CosSatLong = cos(SatLong);    
  189.     LongDiff = SiteLong - SatLong;
  190.     CosBeta = Sil
  191.  reference is in the future.                  */
  192. double GetDay(Year,Month,Day)
  193. double Day;
  194. {
  195.     int TmpYear;
  196.     TmpYear = Year;
  197.     while (TmpYear > RefYear)
  198.    {
  199.    TmpYear -= 1;
  200.    if (TmpYear%4 == 0)
  201.        Day += 366;
  202.     else
  203.        Day += 365;
  204.    }
  205.     while (TmpYear < RefYear)
  206.    {
  207.    if (TmpYear%4 == 0)
  208.        Day -= 366;
  209.     else
  210.        Day -= 365;
  211.    TmpYear += 1;
  212.    }
  213.     if (Year%4 == 0)
  214.         MonthDays[1] += 1;   /* Leap year */
  215.     
  216.      while (Month > 1)
  217.    {
  218.    Month -= 1;
  219.    Day += MonthDays[Month-1];
  220.    }
  221.     if (Year%4 == 0)
  222.         MonthDays[1] -= 1;   /* Leap year */
  223.     return Day;
  224. }
  225.  
  226.